ON PHASE SEGREGATION IN NONLOCAL 
TWO-PARTICLE HARTREE SYSTEMS 



WALTER H. ASCHBACHER AND MARCO SQUASSINA 

Abstract. We prove the phase segregation phenomenon to occur in the ground 
state solutions of an interacting system of two self-coupled repulsive Hartree equa- 
tions for large nonlinear and nonlocal interactions. A self-consistent numerical in- 
vestigation visualizes the approach to this segregated regime. 



1. Introduction 

In this paper, we study the phase segregation phenomenon in the ground states of 
the eigenvalue system consisting of two interacting repulsive Hartree equations whose 
interaction, as well as the respective self-couplings, are nonlinear and nonlocal, 
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Here, the external potentials V\ and V 2 are assumed to be nonnegative and confining, 
whereas the interaction potential V is, for example, of Coulomb type. Moreover, the 
system is purely repulsive, i.e. the self-coupling constants and the interaction strength 
are nonnegative, 

01,02 >0, K>0. 

For fixed i?i,02 > and fixed Ni,N 2 > 0, the phase segregation phenomenon in 
the ground state (0i,02) with ground state energy (^i,// 2 ) of the system (1.1) is 
characterized by the decay to zero of the Coulomb energy functional 

(1.2) D(0i,0 2 )=/ / |0i(x)| 2 \/(x-y)|0 2 (y)| 2 dxdy 
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in the regime of large interaction strength k, i.e. by 

D(0 l5 02 ) — f° r ^ — ^ oo. 

This study is not only of independent mathematical interest but it can also be mo- 
tivated by various physical applications like, e.g., electromagnetic waves in a Kerr 
medium in nonlinear optics, surface gravity waves in hydrodynamics, and ground 
states in Bose-Einstein condensed bosonic quantum mechanical many-body systems 
(see also [1]). The latter domain has been a subject of great interest since many years, 
both on the experimental and the theoretical side, starting off from a series of suc- 
cessful experimental realizations of Bose-Einstein condensation for atomic gases, first 
achieved in 1995 for a single condensate (see e.g. [2]), then, in 1997, for a mixture of 
two interacting atomic species with equal masses (see e.g. [12]), and, finally, in 2003 for 
triplet species states (see e.g. [16]). On the theoretical side, the standard scenario of 
two interacting Bose-Einstein condensates for a very dilute system of repulsive bosons 
uses the description based on a system of two coupled Gross-Pitaevskii equations (see 
e.g. [9, 10, 13, 14]). These equations are formally contained in (1.1) for the case of the 
zero range interaction potential V — 8, i.e. in the case of local nonlinearities. 1 For a 
complete survey paper, we also refer the reader to [7] and references therein. One may 
argue that, for higher density regimes, it is sensible to capture more of the boson-boson 
interaction by allowing for its nonlocal and, hence, less coarse grained resolution by 
use of a potential V ^ 5 (see e.g. [3]). The phase segregation phenomenon has been 
studied e.g. in [13, 15, 17] for Gross-Pitaevskii equations, and it has been given a gen- 
eral variational framework in [6]. Recently, the second author, jointly with M. Caliari, 
has investigated both numerically and analytically the behavior of ground state solu- 
tions 2 highlighting their location 3 and the occurring phase segregation phenomena in 
the highly interacting regime (see [5]). 4 

In the present paper, we extend the analysis to the nonlocal system (1.1) and give a 
proof of the phase segregation phenomenon in the variational calculus setup. Moreover, 
in contradistinction to [5], we adopt a classical self-consistent numerical approach to 
the solution of the ground state of (1.1) in order to compute the phase segregated 
states and to monitor the decay of the Coulomb interaction (1.2). 

As we aim at keeping the paper self-contained and easily readable also for those 
readers who are more acquainted with the physical or the numerical side, we will 
provide rather detailed mathematical arguments throughout the paper. 

2. Strong interaction and phase segregation 

Throughout this section we shall denote by C a generic positive constant which can 
vary from line to line inside the proofs. 

1 For the case of a single condensate, the stationary and dynamical Gross-Pitavskii equation has 
been rigorously derived from the many-body bosonic Schrodingcr equation in the weak coupling limit, 
see e.g. [11] and [8], respectively. 

2 And, in some particular cases, also excited state solutions. 

3 With respect to the off-centering of the confining potentials Vi . 

4 For a complete numerical study of ground states for vector like nonlinear Schrodingcr systems 
with cubic coupling, we also refer to [4]. 
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2.1. Functional setting. As described in the introduction, we are interested in the 
case of nonnegative confining external potentials V\ and Vi- More precisely, we make 
the following assumption. 

Assumption 2.1. The external potentials V\ are nonnegative, continuous, and con- 
fining, i.e. for i = 1,2, we have VJ G C(R N ,M£) with 

lim Vi(x) = oo. 

x|— >oo 

The functional setting we want to apply makes use of the following Hilbert space. 

Definition 2.2. Let the external potentials V{ satisfy Assumption 2.1 and letTi be the 
Hilbert subspace of H 1 ^) x H 1 ^) defined by 

(2.1) H = |(0i,0 2 ) e H\R N ) x H\R N ) : / Vi(x) |0,(x)| 2 dx < oo, i = 1,2), 
where the scalar product of = (0i, fa) G Ti. with ip = (ipi, ^2) Eli. is given fey 

(2.2) (M) H = y2( [ V0,(x) • V^(x)dx+ / ^(x)0^y^(x)dx). 

This functional setting is the natural framework for the study of bound states of sys- 
tems (1.1) in external potentials as it allows (together with Lemma 2.5) the associated 
energy functional (see (2.9)) to be well-defined and finite. 

Lemma 2.3. Under Assumption 2.1, for any 2 < N < 5, the embedding 

is compact, H, being the Hilbert space (2.1) equipped with the norm (2.2). 

Proof. Let (0?, 0^) be a bounded sequence in 7Y, say 11(0?, 02)||h ^ C f° r an h E N. 
Up to a subsequence, it converges weakly in Tt to some (0i,02) G 7Y. Moreover, by 
the Rellich-Kondrachov compactness theorem, up to a further subsequence, (0? , 0^) 
converges strongly to (0i,02) in L 2 (Br) © L 2 (Br) for any i? > 0, where Br denotes 
the open ball in R N of radius R, centered at the origin. Let now M > be an arbitrary 
number. Then, by Assumption 2.1, there exists an R > such that Vt(x) > M for all 
x G \ 5r and any i = 1,2. Hence, we can write 

\tf(x)-<j )t (x)\ 2 dx = [ |^(x)-0 i (a;)| 2 dx+ f \$ (x) - 0,(x)| 2 dx 



M JR N \B R 

f AC 2 
< J b I0f(*)-0,(*)| 2 d* + 1 ^-. 

Let now £ > be given and choose an M > such that 4C 2 /M < £/2. Then, as the 
corresponding radius Rq > is fixed, take h > 1 such that |0f (x) — 0i(x)| 2 dx < 
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e/2 for any h > h . This yields — 4>i\\l 2 - > for h — > oo. Moreover, by the 
Gagliardo-Nirenberg inequality and the boundedness in TC, we have 

Hi - ^ll 4 L1 - < cil^i - M% N M - MbT 2 < CHI - 4n\\% N , 

from which it follows that 

lim - 4n\\ jfa =0. 

This completes the proof. □ 

The interaction between the components fa and fa is described by the following 
Coulomb energy functional which is well-known from classical Hartree theory. 

Definition 2.4. The Coulomb energy functional 5 

D: H\R N ) x H\R N ) -> R 

is defined by 

(2.3) B(fa,fa)= [ [ \fa(x)\ 2 V(x-y)\fa(y)\ 2 dxdy, 



where the interaction potential V is the Coulomb potential in W N for N > 3, 
(2.4) V{x) 



•r 



7V-2 ' 



Due to the following lemma, for any 3 < iV < 6, the Coulomb energy functional 
with potential (2.4) is well-defined. 

Lemma 2.5. Let 3 < N < 6 and let fa e H 1 ^) with \\fa\\ 2 L 2 = N t > for i = 1, 2. 
Then, there exists a constant C such that 



2=2 H=l 
iWh 1 1 1 ^2 1 1^1 ■ 



1 ,fa)<C{N 1 N 2 )^ 
Proof. Due to Schwarz' inequality, we have 

(2.5) B(fa, fa) 2 < B>(fa, fa) B{fa, fa). 

Hence, by the Hardy-Littlewood-Sobolev inequality (for N > 3) and the Gagliardo- 
Nirenberg inequality (for 2 < N < 6), we get 

6-N 

®(fa,fa) < c\\fat L ^ < cii^n^ n^nSr 2 = cN t 2 n^Er 2 , 

which yields the assertion. □ 

Remark 2.6. If N > 2 and the interaction potential is of the form 

(2.6) V\(x) = yrx for some < A < min{4, N}, 
we can again estimate the energy functional 

(2.7) B x (fa,fa):= [ [ \fa{x)\ 2 V x {x-y)\fa{y)\ 2 dxdy 



5 Also called direct term in the Hartree (-Fock) theory. 
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by virtue of the Hardy-Littlewood-Sobolev inequality (for any < A < N) and the 
Gagliardo-Nirenberg inequality (for any < A < 4) as 

BxitiM < C\\M 2 ^Wfof *n < C(iV 1 iV 2 ) i i A ||0 1 ||^ / 1 2 ||0 2 ||^ 2 . 
In particular, if A r = 2 and A = 1, we have 

If A = N — 2 with N > 3 then V>n-2 = © and one recovers the estimate of the previous 
Lemma 2.5. 

2.2. Existence of a minimizer. Let us consider the following two component Hartree 
eigenvalue system in WL N for 3 < N < 6 with Coulomb interaction V from (2.4) and 

N u N 2 >0, 

+ V 1 {x)(j )l + 0! (V * |0i| 2 ) 0! + K (V * |0 2 | 2 ) 01 = ^101, 



(2.8) 



+ V 2 (x)(j)2 + 2 (V * |0 2 | 2 ) 02 + K (V * |0l| 2 ) 2 = yU 2 



2- 



^2\\ 2 L 2=N 2 . 



Since we are interested in the phase segregation phenomenon in the case of a purely 
repulsive Hartree system, we make the following assumption. 

Assumption 2.7. The self-coupling constants 0i,02 and the interaction strength k 
are nonnegative, 

01,02 > 0, K>0. 

Remark 2.8. In the case of coupled Bose-Einstein condensates discussed in the in- 
troduction (where the nonlinearities are local, i.e. V — 5), the self-coupling constants 
0i , 02 as well as the interaction strengths are explicitly related to the scattering lengths 
and the masses of the atomic species in the condensates (see e.g. [9]). 

In order to study the nonlinear ground states of the Hartree system (2.8), we make 
use of the following energy functional. 6 

Definition 2.9. The Hartree energy functional S K : 7i — > [0, oo) defined by 
(2.9) £ K (0i, 2 ) = £oc(0i, 2 ) + «D(0i, 2 ), 

where the decoupled energy functional 8^ : H, — > [0, oo) consists of the sum of the two 
single particle energies Si : 7i — > [0, oo), 



2 



(2.10) £oo(&,&) = 2^£(&), 

i=l 

(2.11) |V0,(x)| 2 dx+ / V l (x)\<f )l (x)\ 2 dx + ^ 

./TUN JmN 2 



6 From here on, since t?i,i?2 > and Ni,N 2 > arc fixed, we display the dependence of the energy 
functionals and the ground state energies on the interaction strength k only. 
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In view of Lemma 2.5, the functional 8 K is well-defined for 3 < N < 6. Moreover, it 
is readily seen that S K is a C l smooth functional and that its critical points constrained 
to the set {(0i,0 2 ) E ||0i|li 2 = ^ f° r * = 1 5 2} are weak solutions of (2.8). 

Remark 2.10. The = corresponds to a noninteracting Hartree system (2.8) 

consisting of two independent Hartree equations, each describing a repulsive single 
particle self-coupling. 

Definition 2.11. The ground state energy E K > of the Hartree functional (2.9) at 
interaction strength n G [0, oo) is defined by 

(2.12) E K = inf 

(<t>i,ct>2)es 

where the infimum is taken over the set 

(2.13) S = {(0i,0 2 ) EH: ||&||£ a = JV 4 /or i = 1,2}. 
Moreover, the segregated ground state energy E^ > is defined by 

E^ = inf £00 (01, 02), 

where now the infimum is taken over the set 

Soo = {(0i,0 2 ) e 5: D(0i,0 2 ) =0}. 

Let us now prove that the Hartree functional (2.9) admits a real and positive mini- 
mizer for any positive interaction strength k. 

Proposition 2.12. Let k G (0, oo) and 3 < N < 6. Then, there exists a positive 
minimizer (0* , 0£) G S of the Hartree functional (2.9) with ground state energy E K 
given in (2.12). 

Proof. In order to prove the assertion, we make use of the direct method in the calculus 
of variations. Hence, we verify the three standard assumptions implying the existence 
of a minimizer. First, since by Lemma 2.3, the normed space 7i from (2.1) with the 
norm from (2.2) is compactly embedded in L 2 (R. N ) © L 2 (R. N ), it follows that the set 
S from (2.13) is weakly closed in 7i. Second, since 

(2.14) ||(01,0 2 )||^ < ^(01,02) 

2 9 

= II (01, fo)\\ 2 H + Ysi D (&> + K 

1=1 

the set {(0i,02) G S : £ K (0i,02) < a} is a bounded nonempty subset of S for any 
positive number a. Third, we have to show that the functional S K is weakly lower 
semi continuous on S. For this purpose, consider a sequence of elements (0^, 2 ) G S 
which converges for h — > oo weakly in 7i to some (0i, 2 ) G S. Since, for any i, j = 1, 2, 

ItfWmm 2 l^*)! 2 !^)! 2 



\ x ~ y\ N 2 \x — y\ N 2 



for a.e. (x, y) G 



q2N 



Fatou's Lemma implies 

(2.15) D(0 1; 2 ) < lim inf D(<#,$), ©(0 4 ,0>) < liminf B(# , 

h—*oo h—>oo 
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Therefore, due to (2.14), the fact that the norm || ■ \\ n from (2.2) on 7i is weakly lower 
semi continuous on S, and (2.15), the Hartree functional £ K is indeed weakly lower 
semi continuous on S. Hence, the three assumptions are verified and the existence of 
a minimizer is proven. Moreover, due to the convexity inequality for gradients, 



|V|0 4 |(x)| 2 dx< / \Vfa(x)\ 2 dx, 
Jr n 



the Hartree functional S K satisfies the following inequality for any (fa, fa) £ «-> ; 

^ K (\fa\, \<h\) < £ K (fa,fa)- 

Consequently, with no loss of generality, we can assume that any minimizer of £ R is 
positive. □ 

Remark 2.13. Note that, for 3 < iV < 5, the Coulomb energy functional D is not 
only weakly lower semicontinuous as given in (2.15), but even weakly continuous over 
H, i.e. for any sequence of elements (fal, 2 ) G S which converges for h — > oo weakly 
in H to some (fa, fa) G S, we have 

(2.16) lim D(0?,$) = B(fa,fa), lim B($,$) = B(fa,fa). 

h— >oo n~>oo 



In order to prove this claim, we make use of Lemma 2.3, which states that the embed- 

4jV 

ding H L ( 
that, for i — 1,2, 



ding TC <^-> L^+a (M^) © (R. N ) is compact. Hence, up to a subsequence, it follows 



(2-17) lim Il0f-0.ll 4^ =0. 

Using (2.17), we want to show that D(0f , 0f ) — > D(0 i; fa) as h — > oo. To this end, we 
use that the Coulomb potential V from (2.4) is even and write 

|D(0f , #) - D(0,, 0,)| < D(| |0f | 2 - I0.H 1 / 2 , (|0f | 2 + |0,| 2 ) 1/2 ). 

By inequality (2.5), the Hardy-Littlewood-Sobolev inequality, and Holder's inequality, 
it follows that there exist a constant C with 

|D(#,#) -D(<fc*)| 2 < C|| ||^| 2 - |*|T /2 |I* ",11 (I* 2 + l<tf) 1/2 ll 4 llif f, 

<C|*f-A||' i]ft . 

This implies, via (2.17), the desired convergence oi3(fa,fa) to B(fa,fa). 7 Hence, 
it follows that all the terms in £ K containing the Coulomb energy functional D are 
weakly continuous on S. Notice that, as a consequence of (2.16), the weak lower 
semi continuity of S K over S also holds in the case of attractive self-coupling or attractive 
interaction, i.e. for i?i,i?2 < or k < 0. In fact, this case amounts to the replacement 
of some (or all) D terms (with positive coupling) in the Hartree functional by — D. 



The convergence D(0i , <t>2 ) — * ID>((^»i , 4>2) as h — > oo can be proved in a similar fashion. 
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2.3. Phase segregation. As pointed out in the introduction, we are interested in the 
situation where the values of the self-coupling constants $i,$ 2 (and Ni, N 2 ) remain 
fixed whereas the interaction strength k becomes very large. 

Definition 2.14. A sequence of minimizers (<fi1,4> 2 ) G S of the Hartree energy func- 
tional S K from (2.9) is said to be phase segregating if 

B(#,$) = o(l) for 00. 

Remark 2.15. If the phase segregating sequence (0J, <p 2 ) is convergent in 7i, then the 
limiting configuration (0^°, 0£°) satisfies D(0~, 0^) = 0. 

Let us now state our main assertion. 

Theorem 2.16. Let 3 < N < Q and let D be the Coulomb energy functional from 
(2.3). Then, for k G (0, 00), any sequence of minimizers (0*, 0£) G S of the Hartree 
energy functional S K from (2.9) is phase segregating for k — > 00, and 

D(0^,0^) = o(^ 1 ). 

In addition, such a sequence converges in the Ti, norm to a minimizer (0i°,02°) G <Soo 
of the decoupled functional Soo from (2.10) and (2.11). 



Corollary 2.17. Under the assumptions of Theorem 2.16, the limiting configuration 
satisfies the following set of uncoupled variational inequalities, 

(2.18) - A0°° + +$i{v* i0r 0r < & 

where N if i°° = %(</>?) + % D(0°°, <PT) and i = 1,2. 



Remark 2.18. Although we stated Theorem 2.16 for minimizers of the Hartree func- 
tional in M N with 3 < A < 6 containing the Coulomb energy functional D, it also holds 
for minimizers of the Hartree functional in M. N with < A < min{4, A} containing 
instead Da from (2.7). This corresponds to the system 



(2.19) 



-A0 1 + V 1 (x)0 1 +^ 1 (V x * 
-A<f) 2 + V 2 (x)<f) 2 + $2(Vx* 



K{V X *\<P 2 \ 2 ) 



/i!01, 
At202, 



110 

110: 



2IIL2 



Ni, 
N 2 , 



where V\ stems from (2.6). In particular, in view of the numerical setup, the case 
A = 2 and A = 1 is covered. 



Proof. Consider a sequence of minimizers (0J , <fi 2 ) G S for k — > 00 whose existence 
is assured by Proposition 2.12. Note first that, in the light of Definition 2.11, the 
sequence of corresponding ground state energies {E K ) is uniformly bounded because 

(2.20) E K = inf £«(0i,0 2 )< inf ^(0 1; 2 )= inf 8 00 {<f>u<h)=E 00 . 

{<t>l,<t>2)&S Wl^2j6Soo (01,<fe)6Soo 
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In particular, due to (2.14) and the definition of a minimizer, the sequence (0£, (f>% ) is 
uniformly bounded in Ti with respect to k, 

\m,^)\\h<£MA2) = e k <e 00 . 

Hence, since TL is weakly sequentially compact, there exists a pair (4>f, 2 x> ) G H and a 
subsequence of ((f)*, 0£), again denoted by (0£, 0£) which, for k — > oo, converges weakly 
in H to ((f)f,(f)f). Next, we want to show that (^Af) e ^oo- Since G 5 

and the embedding L 2 (M Ar ) © L 2 (M 7V ) is compact, we have, for i — 1,2, 



LOO II 2 



NVi IlL 2 - -^i- 

Hence, (0i°,02°) G S. Moreover, again due (2.14) and (2.20), we have 

(2.21) ^©(01, 02) < < Eoo, 

which implies that the sequence (0*, 0£) is phase segregating for k — > oo, 

0(0?, 0*) =0(O- 

Also, since we know that the Coulomb energy D is weakly continuous on S, we get 

o(0r\0-) = o, 

and, therefore, (0i°,02°) G <Soo- In order to prove that (0^,0^°) is a minimizer of 
and that (0J, 0^) converges strongly in H to (0£°, 0^), we next show that D(0i, 05) = 
o(k~ 1 ). To this end, consider the sequence kD(0^, 0^) which is bounded due to (2.21), 
and pick a convergent subsequence, denoted by K n D(0" n , 02™) • Then, using that 
(0i° 5 02°) e ^oo, the weak lower semi continuity of the decoupled energy functional 
on S, and (2.20), we get 



(2.22) £oo(0r,0n + lim ^0(0^,0*") < liminf £ Kn ( 



<^oo<^oc(0 1 °,O, 

with the consequence that K n D(0^ n ,02™) = as n — > oo. Therefore, since this holds 
for all convergent subsequences of k 0(0?, 05), we arrive at 

(2.23) D(0i,02) =o(«" 1 ). 

This implies, on one hand, that (0?, 0£) converges strongly in 7i to ((f)f, (f>??) since 
from E k < ^00(0^°, 2 >o ), (2.23), and the weak continuity of D(0f , 0f), we get 

limsup||(0^,02 K )||i< 11(0^02°) 111- 



On the other hand, using (2.22) and (2.23), we see that (0i°,02°) is a minimizer of 
Soq, that is Eoq = £00(01°, (f)^). Finally, we note that, again due to (2.22), we have 
E K — > i?oo as k — > 00. This brings the proof of Theorem 2.16 to an end. □ 

Finally, we prove the assertion of Corollary 2.17. 
Proof. Observe that, by virtue of 

(2-24) = ±{£iW) + y B(#, tf) + kD(#,# 
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the inequality (2.21) and £i((f>f) < E^, we get 

sup /i^ < oo, 

where denotes the nonlinear eigenvalue of the minimizer 0^ as the weak nonlinear 
ground state in the corresponding nonlinear eigenvalue system (2.8). Then, up to a 
subsequence, \i1 — > as k — > oo. Testing the equations of (2.8) with arbitrary 
nonnegative functions r\ of compact support, we get, recalling that <pi > 0, 

/ V4>i(x) - Vrj(x) dx + / Vi(x) (f>i(x)r](x)dx 
Jr n Jr n 

+*[ f Mfm^w dxdy 

Jr n Jr n \ x — y\ 

< / 77(x) dx. 

Hence, letting «; — > oo, it turns out that 0°° satisfies the variational inequality (2.18). 
Finally, the strong convergence and (2.24) yields, for % — 1,2, 

This ends the proof of Corollary 2.17. □ 



3. Numerical approach 

3.1. Galerkin approximation of the nonlinear eigenvalue system. In order to 
carry out the numerical simulation, we treat the Hartree system in the plane from 
(2.19) with N = 2 and A = 1 in the framework of the following finite element approx- 
imation. As physical subdomain of M 2 , we choose the open square 

(3.1) n=(0,D) x2 

with D > whose closure is the union of the (m — l) 2 congruent closed subsquares 
generated by dividing each side of Q equidistantly into m— 1 intervals. Let us denote by 
M = (m — 2) 2 the total number of interior vertices of this lattice and by h — D/(m— 1) 
the lattice spacing. 8 Moreover, let us choose the Galerkin space Sh to be spanned by 
the bilinear Lagrange rectangle finite elements tpj G C(Q). 9 
Hence, with this choice, we have 

S h C C(H) n^fi) 

8 As bijection from the one-dimensional to the two-dimensional lattice numbering, we may use the 
mapping r : {0, m — l} x2 — > {0, m 2 — 1} with j = r(mi, m 2 ) := mi + ra-im. 
9 The reference basis function tpa : — > [0, oo) is defined on its support [0, 2h] x2 by 

xy, if (x,y) € [0,/i] x2 , 

{2h-x)y, ii(x,y)e[h,2h]x[0,h], 

(2h-x)(2h-y), if(x,y)e[h,2h]* 2 , 

x{2h-y), ii(x,y)e[0,h]x[h,2h], 

sec Figure 1. The functions (pj arc then defined to be of the form (3.2) having their support translated 
by (mi/i, m^h) with m 1; m 2 = 0, m — 3. 



(3.2) 



<A) (z, y) 
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Figure 1. (po(x,y) on its support [0,2/i] with maximum at vertex (h,h) 



and dimS'/j = (m — 2) 2 . The Hartree system (2.19) in its weak finite element approxi- 
mation form reads, for all (p e 5^, 10 

f (Vy?, V0i) + (<p, V x <f>i) + (<p, (V * [VM 2 
(V(p, V0 2 ) + (cp, V 2 <f) 2 ) + (cp, [V * [tf 2 |0 2 | 2 

\\Mv = Ni, 
U2\\h=N 2 . 



(3.3) 



«|02| 2 ])0l) = A*1,0 (Vi^l) 
«#l| 2 ])02) = A*2,o(y»02), 



If we expand a G 5/, with expansion coefficients z c 
finite element basis {^j}^ of Sh, 



[z a ,l, .-,Z a ,M) e CM w - T - t - tne 



(3.4) 



and if we plug this expansion together with y? = cpi, i 
following coupled matrix system on C M , 



M, into (3.3), we get the 



(3.5) 



Ql[zi,Z 2 ]zx — /Xo.l^l) 
Q2[Zl, Z 2 ]z 2 = fiO, 2 Z 2 , 

iV 2 . 

Here, the matrix- valued mappings Qi,Q 2 '■ C M x C A/ — > C are defined by 

Q 1 [z 1 ,z 2 ] := A-^S + n + ^G^il+KG^]), 
Q 2 [*i,2 2 ] := A- 1 (B + Y 2 + ^ 2 G[z 2 ]+kG[z 1 \), 



\v I 2 
1^1 1 2 

|y I 2 
I P2|2 



10 In this section. (■, ■) and || • || stand for the L 2 (f2)-scalar product and L 2 (S7)-norm, respectively. 
Moreover, the convolution on the finite domain O is defined, for any (x,y) € f2, by (V * <f>)(x,y) := 
J n V{x-x',y-y')^{x',y')dx'Ay'. 



11 



(z, z)i denotes the L -norm on 



where 



Sjli and ( z : w h 



i,Aw) 



are the Euclidean and the L 2 -scalar product on C M , respectively. For </> a from (3.4), we have ||^a|| 2 = 



6 a 12 



iV Q for a = 1,2. 
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and A e C MxM is the mass matrix, B G C MxM the stiffness matrix, and Y a e C MxM 
the matrices generated by the external potentials V a , 

(3.6) Aij := (<Pi,(pj), Bij := (VipiiVtpj), (Ya)ij ■= (fu V a (pj). 

Moreover, the matrix- valued mapping G : C — > C is defined oiaw = (w\, wm) € 
C M by 

(A/ \ M 

ipi,g[y^W k ip k ] ipj J = ^ W k W l V iklj, 
k=l / k,l=l 

where the function g and the Hartree convolution term V^ij are defined by 

(3.7) g[<j>] :=V* |</>| 2 , V m ■= [tfu V * 

Remark 3.1. We avoid the inversion of the mass matrix A and simplify the evaluation 
of the double integral in the Hartree convolution term (3.7) by approximating the 
integrals over Q by the standard mass lumping quadrature procedure. 

In order to simplify the eigenvalue system (3.5) with the help of Remark 3.1, let us 
introduce the mappings H 1 ,H 2 : C M x C M -> C MxM defined by 

#i[*i,z 2 ] := ^(fi + Vl +^ 1 diag(G [^i]) + «diag(G h])), 

H 2 [z u z 2 ] := — ( J B + y 2 + ^ 2 diag(G h]) + Kdiag(G , oki])), 

where diag : C M — > C is defined to be the matrix- valued mapping on w = 
(wi, ...,Wm) £ C M defined by diag(w)jj := SijWj for all i,j = 1, M, and Go : C A/ — > 
C M is defined by 

M 
3=1 

where r is the grid numbering bijection from footnote 8. Hence, Remark 3.1 amounts 
to the replacement G[z] i— > diag(Go[z]) and we get approximated Hartree system 

, Z 2 \Z\ — /io.1^1) 
Z 2 \ Z 2 = ^0,2^2, 

\zi\i = Ni, 
[\z 2 \ 2 2 = N 2 . 



(3.8) 



3.2. Algorithms. In order to solve the nonlinear coupled eigenvalue system (3.8), we 
make use of the method of successive substitution 12 whose fixed-point map is con- 
structed with the help of the power method used for the solution of the corresponding 
linearized problem. In the following, we briefly describe the basic ideas of these algo- 
rithms. 



Also called nonlinear Richardson iteration or Picard iteration. 
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Method of successive substitution (MSS) 
Let M. be the compact set 

M := {[zi, z 2 \ E C AI x C M | \ Zl \l = N u \z 2 \\ = N 2 }. 

The MSS is an iterative method of the form 

(3.9) [zf +l \zt +l) ]=F[zf\zt\ 

where the fixed point map F : M. — > M. is constructed as follows. Given an approx- 
imate nonlinear system ground state [z^, G M. at iteration level n G N, the 
approximate nonlinear system ground state [^ n+1 , z 2 +1 ] G M. at iteration level n+ 1 
is defined to be the linear system ground state of the linearized eigenvalue system 

(") Ji)l_(n+1) _ in+1) (n+1) 



(3.10) 



o.i z i i 
o- \Jp) _(n)i_(n+l) _ (n+1) (n+1) 
^2^1 , ^2 J z 2 _ e 0,2 z 2 J 

\ z l \2 ~ iV l' 
I (n+1) 1 2 _ at 
\ Z 2 \2 — iV 2- 



in) »i 



Remark 3.2. Here and in the following, we make the assumption that H a [z\ ,z 2 
has a unique linear ground state. E.g., using perturbation theory in the regime of 
small nonlinear couplings, this holds as soon as the linear operator if Q [0,0] has a 
nondegenerate ground state energy. 

Remark 3.3. We can write the fixed point map F a [zi l \ z 2 1 ^} from (3.9) with the help 
of the linear ground state projection 

P a [zi,z 2 ] = ~ I (H^z^-C)' 1 dC, 



2vri 



r Q [21,22] 



where r a [zi, z 2 ] is a path which encircles the linear ground state energy of H a [z±, z 2 ] 
in the positive direction and no other point of the spectrum of H a [z±, z 2 \. The map F 
can now be written as 

(3.11) F a [z^A n) ] = V^a „ „ m ■ 

V 1 OLi'iJVa („) („) ( n ) 

Since H a [-, ■} is Lipschitz continuous on Ai, the map F has a not necessarily unique 
fixed point due to Schauder's fixed point theorem. 

The system (3.10) being not only linearized but also decoupled, we can solve the 
two linear eigenvalue problems separately. In order to approximately determine the 
ground states of the linear eigenvalue problems, we make use of the power method 
which works as follows. 

Power method (PM) 
The PM computes the eigenvector of H^z^, z^] whose eigenvalue has largest modu- 
lus amongst all the eigenvalues whose eigenvectors appear in the eigenvector expansion 
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of the starting approximation. To access the ground state of H a [zi , ], we apply 
the following spectral shift 13 

:= \H a [zt\g% + l. 
Moreover, we define the shifted operator by 

H a [z[ n \4 n) ] :=H a [zt\zf ] ]-s^. 
Now, the p-th iterate of the PM iteration is defined by 14 

H \z {n) z (n) fz {n) 

/oi 9 \ Jn+l),p ._ n o\ z \ i Z 2 J Za 

1 ^ ' \tt r (n) (n)f (n),- 

| -"a |*1 5 z 2 \ Za I 

Remark 3.4. Since H a [zi , z^ ] is real symmetric, the spectral theorem implies the 
existence of an orthonormal basis of C M of eigenvectors {wa^}^ 1 ^ 1 of H a [z[ n \ 2^ ]- 15 
Moreover, since the spectral radius of H a [z^{ \ z^\ is smaller than s£\ we have for all 
eigenvalues of the shifted operator e a ^ € spec(H a [z[ n \ z^]) — s« that 

-2s^ ] < e Q , < e Q ,i < ... < e a>M -i < 0. 

Let us expand z^ 1 w.r.t. the orthonormal basis {w^k}^^ 1 as 

M-l 

Z a ^ ^ , ^,a,k^a,ky 
k=0 

use that H a [z[ , z\ ]u> a ,& = e a ,kW a ^, and divide the numerator and the denominator 
in (3.12) by |e a ,o| p - Moreover, let us assume that £ a) o ^ 0. Then, in the large p limit, 
(— l) p z^ +l ' >,p converges to a multiple of the ground state of H a [z[ n \ z% ], 

(3.13) 4»+i)* = {-If w 0>a + o(l). 

|s0,a| 

Remark 3.5. Using formula (3.13), the fixed point map F from (3.9), (3.11) can also 
be written as 

1 2a 1 2 

Stopping criteria 

Both for the inner PM iteration and the outer MSS iteration, we use a relative error 
stopping criterion in the numerical computation. For the PM iteration, let us define 
the energy 

(3.14) e to +1) ' P := (4 n+1) ' P , H a [z[ n \ 4 n) ]4 n+1) ' P >- 



13 For A = [ciij] e C MxM , we define the I 1 - matrix norm by \A\ 1 := X^-=i l a u"l- 
14 |z| := (z, z) 1 / 2 denotes the Euclidean norm of z € C M . 

15 We suppress the superscript n in the eigenvectors, eigenvalues, and in the expansion coefficients, 
since the PM iteration acts at a fixed n. Moreover, the numbering starts at being the index of the 
ground state. 



ON PHASE SEGREGATION IN NONLOCAL TWO-PARTICLE HARTREE SYSTEMS 15 

Then, for suitably chosen accuracy tolerance <5pm > 0, we stop the PM iteration for 
each component a = 1 , 2 as soon as 

\{TT uM ~( n )i ;(™+l),P\ Jn+l),P\ 
/„ \\ I1 a[ Z l ,Z 2 J t aQ ) Z a | 

14,0 I 

Remark 3.6. Note that the quotient (3.15) does not depend on the shift s« , since 
the iterates Za 1 p are normalized w.r.t. the Euclidean norm on C M . 

Remark 3.7. Clearly, the stopping criterion (3.15) is satisfied for any eigenvector of 
H a [z^,Z2 ]. But as soon as £ Qi0 7^ 0, e.g. due to finite precision arithmetic, the PM 
iterate z^ +1 ^' p converges to a multiple of the ground state w a $. But note that the 
chosen accuracy may be reached before a nonvanishing ^ Q) o is generated. 

For the MSS iteration, we implement a similar stopping criterion. To this end, we 
define the approximate nonlinear ground state energies as 

„ ■- 1 /,(n+l) tr rjn+l) _(»+l)l-(n+l)\ 
^0,0 - — ir \ z o ;-"a|*l ) z 2 J z a /2; 

where, compared to (3.14), the Hartree energy H a depends on iteration level n + 1 
instead of level n. We stop the MSS iteration as soon as 



iH a [zt +i \zt +i) ]-^: l) )^ +1) \ 



2 -<5 



MSS; 



where #mss > is some suitably chosen accuracy tolerance. 

3.3. Phase segregation. As it has been defined above in Definition 2.14, a sequence 
of nonlinear ground state solutions (0^, is phase segregating if its Coulomb energy 
vanishes in the limit of large interaction strength k, i.e. 

(3.16) D(^,^) = (^,(\/*|02l 2 )^)^O for/^oo. 
Plugging the expansions (3.4) into (3.16), we get 

M M 
i=l j=l 

Hence, making use of Remark 3.1, we define the approximated Coulomb energy Do : 
C M x C M — > R by 

(3.17) B [zi, z 2 ] := (zi, diag(G [«2])^i). 

Below, we will use this approximation in the numerical computation of the Coulomb 
energy. 
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4. Figures 



The numerical computations leading to the following figures visualize the qualita- 
tive picture of the approach to the segregated regime. First, we exhibit the densities 
of the wave functions <pi and 02 for increasing values of the interaction strength k 
approaching the segregated regime. Second, we report on the decay of the Coulomb 
energy (3.17). 

We choose the external potentials V a for a = 1, 2 to be isotropic harmonic potentials, 
(4.1) V a (x, y) = c a ((x - a a f + (y - b a ) 2 ) , 

and the interaction potential V to be a regularized Yukawa potential, 



(4.2) 



V(x,y) 



yjx 2 + y 2 + 7 

Remark 4.1. The potential (4.2) being the regularized three-dimensional Yukawa 
potential, it may be argued that we consider a physically three-dimensional system 
constrained to a two-dimensional submanifold of the three-dimensional configuration 
space. 

The specification of the parameters used in the simulations below is summarized in 
the following table (cf. (3.1), (3.3), (4.1), and (4.2)). 16 





N 2 


ai 


h 


Cl 


a 2 


b 2 


c 2 


A 






r 


7 


1 


1 


D/2 


D/2 


10 5 


D/2 


D/2 


10 3 








cf. below 


10 2 


KT 1 



Remark 4.2. All the qualitative features of the following simulations have been tested 
for stability in different physical and numerical parameter ranges. 

4.1. n = 0. For the interaction strength k = 0, the system is uncoupled and linear, 
and we find the ground state wave functions of the harmonic oscillator. The supports 
are fully overlapping, see Figure 2. 17 

4.2. n = 0.5. The wave functions 0i and 2 start to feel their respective repulsion. 
The support of 0i is retracting whereas the one of 02 gets pushed outwards. The 
supports are still heavily overlapping, see Figure 3. 



The code is part of our Hartree package written in C++. 
17 A11 the fi gures have been produced with the help of gnuplot. 
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FIGURE 3. The interaction strength is k = 0.5. 
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4.3. k = 10. In the regime of large interaction strength k, the segregation phenomenon 
occurs: the supports get more and more disjoint, see Figure 4. 




Figure 4. The interaction strength is k = 10. 

Remark 4.3. Up to the shape of the support of 0j, there is no qualitative change in 
the picture if the two harmonic potentials are slightly dislocated with respect to each 
other. 

4.4. Coulomb energy. Finally, we monitor the decay of the Coulomb energy from 
formula (3.17), see Figure 5. 




2 4 6 8 10 



Figure 5. The decay of Dofci , z%\ as a function of k. 
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